/*
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
#if !defined _OS_American_Pricer_
#define _OS_American_Pricer_

#include <src\algorithm\LinearAlgebra/TridiagonalSystem.h>
#include <src\algorithm\Evolution_Algorithm\all.h>
#include <src/math/CauchyProblems/CauchyProblem_traits.h>
#include <src/math/BoundaryCondition/BoundaryCondition.h>
#include <src/math/Stochastic/Stochastic_Tree.h>
#include <src/math/Stochastic/StochasticProcesses.h>
#include <src/miscellaneous/Function_Utilities.h>



/** @ingroup American_Pricer
*   @brief an Euler scheme to solve the free boundary problem @link American_Pricer@endlink
*
*/
template <class tree_type, class payoff_type>
	class OS_American_Pricer_Euler : public Euler<
		Delta_t<typename tree_type::underlying_vector_type, 
		tree_type::underlying_vector_type::Dim>
										> 
	{
	public :
		typedef Euler											Base;
		typedef tree_type										tree_type;
		typedef typename tree_type::stochastic_process_type		stochastic_process_type;
		typedef typename tree_type::underlying_vector_type		underlying_vector_type;
		typedef typename tree_type::density_vector_type			density_vector_type;
		typedef payoff_type										payoff_type;
		
		OS_double maturity_;
		OS_double from_time_; 
		OS_double to_time_;
		OS_double time_12_;
		OS_double x0_;
		OS_double compound_; 
		OS_double drift_;
		OS_double volat_;
//		OS_double strike_;

		OS_American_Pricer_Euler() : steps_(0) {};
		void set_parameters(smart_ptr<tree_type> & tree,smart_ptr<payoff_type> & payoff) {
			tree_ = tree;
			payoff_ = payoff;
			maturity_ = payoff_->get_maturity();
//			strike_   = payoff_->get_strike();
			x0_ = tree_->get_process()->initial_conditions();
		};
		virtual OS_double next_time(time_type to) 
		{
			Base::next_time(to);
			from_time_ = maturity_-from_;
			to_time_ = maturity_-to_;
			time_12_ = (from_time_ + to_time_) /2.;
			stochastic_process_type & process = *tree_->get_process();
			OS_double rate  = process.rate(to_time_,x0_);
			compound_  = std::exp(-rate*delta_t_);
			drift_ = process.OS_Base::drift(to_time_,x0_); // this is (r_1/2-\sigma ^2)
			volat_ = process.OS_Base::volatility(to_time_,x0_);
			return to_;
		};

		virtual OS_double CFL() {
			return maturity_ / (steps_);
		};

		virtual void vstar() {};
		virtual void gstar() {};
		virtual void update(){
				u_ = *get_state_variable();
				steps_ = tree_->get_steps();
		}; 

		virtual void operator ++() {
//			OS_double test = jmmath::moment1(u_);

			underlyings_ = tree_->underlying(to_time_); // these are the sampling points at time from_time_

			underlying_vector_type S(tree_->inverse_cumulative_);
			OS_Array_multiply(S,volat_ * std::sqrt(from_time_));
			OS_Array_plus(S,drift_*from_time_);
			OS_Array_exp(S);
			OS_Array_multiply(S,x0_);

			cumulative_ = OS_cumulative_normal_distribution<OS_double>(drift_*delta_t_,volat_*std::sqrt(delta_t_));
			underlying_vector_type local(steps_);
			local.assign(0.);

			underlying_vector_type probatest(steps_);
			probatest.assign(0.);

			typedef underlying_vector_type::iterator iterator;
			iterator underlyings_begin = underlyings_.begin();
			iterator local_begin = local.begin();
			iterator underlyings_end = underlyings_.end();
			for (;underlyings_begin != underlyings_end; ++ underlyings_begin, ++local_begin)
			{
				OS_double Pi = 0;
				OS_double x_0 = *underlyings_begin;
				iterator u_begin = u_.begin();
				iterator S_begin = S.begin();
				iterator S_end = S.end();

				iterator proba_begin = probatest.begin();
				for (;S_begin != S_end ; ++u_begin, ++S_begin, ++proba_begin)
				{
					OS_double cumulative = Cumulative(x_0,*S_begin);
					OS_double proba = cumulative-Pi;
					*local_begin += proba * (*u_begin);
					*proba_begin += proba;
					Pi = cumulative;
				};
			};

//			std::cout << std::endl << probatest << std::endl;
//			std::cout << std::endl << jmmath::norm1(probatest) << std::endl;
			u_ *= compound_;
			Exercise();

//			std::cout << std::endl << (jmmath::moment1(u_) - test)/test << std::endl;
		};
		virtual void Exercise()
		{
				underlying_vector_type underlying(tree_->underlying(to_time_));
				solution_type payoff ( (*payoff_).operator()<underlying_vector_type,solution_type>(underlying,to_time_) );
				OS_Array_max(u_,payoff);
		};

		OS_size steps_;
		smart_ptr<tree_type>					tree_;
        smart_ptr<payoff_type>					payoff_;
		underlying_vector_type underlyings_;
//		underlying_vector_type inverse_cumulative_; // we store S_{\mathcal{N}}(1,y_{i,i+1})
//		density_vector_type normal_densities_;		// we store \mathcal{N}(1,S_{\mathcal{N}}(1,y_{i,i+1})) 
	private :
		OS_cumulative_normal_distribution<> cumulative_;
		inline OS_double Cumulative(OS_double x_0, OS_double x) {
			return cumulative_(std::log(x/x_0));
		};
	};

	template <class tree_type, class payoff_type>
	class OS_American_Pricer_Base_CN : public OS_American_Pricer_Euler<tree_type, payoff_type>
	{
	public :

		typedef OS_American_Pricer_Euler					Base;
		underlying_vector_type						underlying_from_time_;  //underlying at previous time (from_time_)
		underlying_vector_type						underlying_to_time_;	//underlying at next time     (to_time_)
		solution_type								payoff_to_time_;		// payoff values at next time
		solution_type								u_1_;					// solution at previous time
		density_vector_type							rhs_;					// second member of the tridiagonal linear system
		density_vector_type							upper_diagonal_;		// diagonal_ of the tridiagonal linear system
		density_vector_type							diagonal_;				// upper diagonal_ of the tridiagonal linear system

		virtual OS_double next_time(time_type to) {
			Base::next_time(to);
			if (std::abs(from_time_ - maturity_) < THRESHOLD) 
				underlying_from_time_ = tree_->underlying(from_time_);
			else
				underlying_from_time_ = underlying_to_time_;
			underlying_to_time_ = tree_->underlying(to_time_);
			payoff_to_time_ = (*payoff_).operator()<underlying_vector_type,solution_type>(underlying_to_time_,to_time_ );
			return to_;
		};
		virtual void Exercise()
		{
				OS_Array_max(u_,payoff_to_time_);
		};
		virtual void update(){
			Base::update();
			densities_.resize(steps_);
		};
		virtual void compute_densities(OS_double from_time_, OS_double to_time_) 
		{
			underlying_vector_type S_n(underlying_to_time_);
			OS_Array_divide(S_n,compound_);

			densities_ = (underlying_from_time_+S_n);
			densities_ = 4./ ( ( (densities_<<1) - densities_ )*steps_ );

			typedef density_vector_type::iterator density_iterator;
			density_iterator density_begin = densities_.begin();
			density_iterator density_end = densities_.end();
			typedef underlying_vector_type::iterator undelying_iterator;

			undelying_iterator S_n_1_iterator = underlying_from_time_.begin();
			undelying_iterator S_n_iterator = S_n.begin();
 

			OS_double integrator(0.);
			for (;density_begin != density_end;++density_begin,++S_n_1_iterator,++S_n_iterator)
			{
				integrator += ( (*S_n_iterator) - (*S_n_1_iterator) )/delta_t_/steps_;
				(*density_begin) *= integrator;
			};
/// test positivity of density. Note : positivity may be true only under a CFL condition
//			bool chek_positivity = OS_Array_less<OS_double,density_vector_type>(0.,densities_);
//			if (!chek_positivity)
//				OS_DYNAMIC_CHECK( false, " negative densities"  ); //test the moment conservation
///

		};
		virtual density_vector_type second_member(underlying_vector_type & initial_condition) {
			density_vector_type result(initial_condition);
			result = ( (result<<1) - result )*(steps_);   // \delta_y f^{n+1}
			OS_Array_xor(result,densities_);    // a_i \delta_y f^{n+1}
			result =  (result - (result>>1) )*( delta_t_*steps_/4.); // \frac{\tau^n}{2} \delta_{\bar{y}} a_i \delta_y f^{n+1}
			result += initial_condition; // (Id + \frac{\tau^n}{2} \delta_{\bar{y}} a_i \delta_y ) f^{n+1}
			return result;
		};
		virtual void compute_tridiagonal_operator() 
		{
			diagonal_ = densities_;
			upper_diagonal_ = densities_;
			diagonal_ += densities_ >>1;
			OS_double local_constant (steps_*steps_*delta_t_/4.);
			OS_Array_multiply(upper_diagonal_,-local_constant);
			OS_Array_multiply(diagonal_,local_constant);
			OS_Array_plus(diagonal_,1.);
		};


		virtual void operator ++() {

//////////////////test purposes : stores the first moment of the solution at previous time
//			u_1_ = u_;
//			OS_double test = jmmath::moment1(u_); // test
////////////////////// test

			compute_densities(from_time_, to_time_);

// second member is stored in rhs_
			rhs_ = second_member(u_);
			
////////////////////// test the first moment conservation of the right hand side
//			double check_test = std::abs(test - jmmath::moment1(rhs_))/test;
//			if (check_test > 0.00001)
//				OS_DYNAMIC_CHECK( false, " first operator is not conservative"  ); //test the moment conservation
////////////////////// test

// tridiagonal operator
			compute_tridiagonal_operator();
			
////////////////////// test the tridiagonals operators
//			density_vector_type local_test = Apply(u_,diagonal_,upper_diagonal_,upper_diagonal_); //test
//			check_test = std::abs(test - jmmath::moment1(local_test) )/test;
//			if (check_test > 0.00001)
//				OS_DYNAMIC_CHECK( false, " second operator is not conservative"  ); //test the moment conservation
//			local_test = TridiagonalSolver(local_test,diagonal_,upper_diagonal_,upper_diagonal_);
//			check_test = jmmath::norm1(u_-local_test);
//			if (check_test > 0.00001)
//				OS_DYNAMIC_CHECK( false, " problem with the tridiagonal solver / application"  );					//test the trilinear solver
////////////////////// end test

			u_ = TridiagonalSolver(rhs_,diagonal_,upper_diagonal_,upper_diagonal_); // this is the linear part
////////////////////// test the first moment conservation
//			check_test = std::abs(test - jmmath::moment1(u_) ) / test;
//			if (check_test > 0.00001)
//				OS_DYNAMIC_CHECK( false, " first moment is not conservative in CN"  ); //test the moment conservation
////////////////////// test
////////////////////// test the densities computations
//			local_test = second_member(underlying_from_time_);
//			local_test = TridiagonalSolver(local_test,diagonal_,upper_diagonal_,upper_diagonal_); // this is the linear part
//			OS_Array_multiply(local_test,compound_);
//			check_test = jmmath::norm1(local_test - underlying_to_time_) / jmmath::norm1(local_test);
//			if (check_test > 0.00001)
//				OS_DYNAMIC_CHECK( false, " there is a problem computing the transition probabilities"  );
//////////////////////

			OS_Array_multiply(u_,compound_);
///////////////////// test the first moment of the linear part of the solution
//			check_test = jmmath::moment1(u_)/compound_;
//			check_test = ( std::abs(check_test - test)) / test ;
//			if (check_test > 0.00001)
//				OS_DYNAMIC_CHECK( false, " first moment evolution is not good"  );
//////////////////////////////////////////////////////////
			Exercise();

		};
		density_vector_type densities_;
	};

	template <class tree_type, class payoff_type>
	class OS_American_Pricer_CN : public OS_American_Pricer_Base_CN<tree_type, payoff_type>
	{
	public :
		typedef OS_American_Pricer_Base_CN									Base;
		typedef typename underlying_vector_type::iterator					underlying_iterator;
		typedef std::pair<underlying_iterator,OS_size>				iterator_pair_type;
		typedef std::list<iterator_pair_type>						iterator_pair_list_type;
		typedef typename iterator_pair_list_type::iterator			iterator_pair_list_iterator_type;
		iterator_pair_list_type										exercised_list_;

		virtual void Exercise()
		{
			exercised_list_.clear();
			typedef solution_type::iterator payoff_iterator;

			payoff_iterator begin_payoff = payoff_to_time_.begin();
			payoff_iterator end_payoff = payoff_to_time_.end();

			payoff_iterator begin_u = u_.begin();

			exercised_list_.push_front(iterator_pair_type(begin_u,0));

			bool exercised(true);

			for (;begin_payoff != end_payoff; ++begin_payoff, ++begin_u) {
				if (*begin_payoff > *begin_u + THRESHOLD) 
				{
					*begin_u = *begin_payoff;
					exercised = true;
				}
				else {
					if (exercised) {
						exercised_list_.push_front(iterator_pair_type(begin_u,0));
					};
					++(exercised_list_.begin()->second);
					exercised = false;
				};
			};
		};

		virtual void operator ++() {
			Base::operator++();
			for (iterator_pair_list_iterator_type begin = exercised_list_.begin();
				begin != exercised_list_.end();
				++begin) 
			{
					OS_double local_constant (steps_*steps_*delta_t_/4.);

					OS_size left	= begin->first - u_.begin();
					OS_size right	= begin->second;
					OS_size size	= right;
					right += left;
					if ((size > 1) && (size < steps_)) //TODO the case where only one point is exercised
					{
						if (left > 0)
							rhs_[left] += u_[left-1]*(densities_[left-1]*local_constant/compound_);
						if (right < steps_)
							rhs_[right-1] += u_[right]*(densities_[right]*local_constant/compound_);
						TridiagonalSolverIterator<solution_type,underlying_vector_type>
						(
							u_.begin()+left,
							rhs_.begin()+left,
							diagonal_.begin()+left,
							upper_diagonal_.begin()+left,
							upper_diagonal_.begin()+left,
							size
						);
						OS_Array_multiply(u_.begin()+left,u_.begin()+right,compound_);
					};
			};


/////  test the solution
////	condition (f >= P)
/*
			Exercise();
			bool check_bool = exercised_list_.begin()->second == steps_; // means no exercising
			if ( !check_bool)
				OS_DYNAMIC_CHECK( false, " the payoff is greater than the solution"  );
////	condition (f-P)*(\partial_t f - \frac{1}{2}\partial_y \sigma^2\rho^2\partial_yf)=0
			underlying_vector_type test_u(u_ / compound_);
			density_vector_type local_u((u_1_ + test_u)*.5);
			local_u = ( (local_u<<1) - local_u )*(steps_);
			OS_Array_xor(local_u,densities_);
			local_u =  (local_u - (local_u>>1) )*( steps_/2.);
			test_u = (u_1_-test_u)/delta_t_ + local_u;
			std::cout << std::endl << test_u;
			underlying_vector_type verif(u_-payoff_to_time_);
			std::cout << std::endl << verif;
			verif ^= test_u;
			OS_double check_test = jmmath::norm1(verif);
			if (check_test > 0.00001)
				OS_DYNAMIC_CHECK( false, " the condition (f-P)*(partial_t f - frac{1}{2}partial_y sigma^2 rho^2 partial_yf)=0 is not verified"  );
*/
////end test
		};
	};


	template <class tree_type, class Payoff, class one_step_algorithm>
	class Cauchy_American_Pricer : public CauchyProblem_traits<one_step_algorithm,Graphical_Interface>
	{
		typedef Unbounded_Vector<double> underlying_vector_type;
		typedef Dirichlet_Vector<double> density_vector_type;
		typedef typename tree_type::stochastic_process_type		stochastic_process_type;
		public :
			Cauchy_American_Pricer(smart_ptr<tree_type> &tree,smart_ptr<Payoff>& payoff) : payoff_(payoff), tree_(tree)
			{
				get_algorithm().set_parameters(tree_, payoff_);
			};
			virtual void set_dim(OS_size dim) {
				OS_DYNAMIC_CHECK(dim == 1,"Cauchy_StochastiqueOnly Dim 1");
			};
			OS_size get_dim() {return 1;};
			void initial_conditions(OS_size steps) {
				tree_->set_steps(steps);
				initial_conditions();
			};
			void initial_conditions()
			{
				double maturity_ = payoff_->get_maturity();
				underlying_vector_type underlyings(tree_->underlying(maturity_));
				density_vector_type payoff ( (*payoff_).operator()<underlying_vector_type,density_vector_type>(underlyings,maturity_) );
				set_initial_conditions(payoff);
			};
			virtual smart_ptr<interface_type> Create_Instance() {
				return smart_ptr<interface_type>(new Cauchy_American_Pricer(tree_, payoff_));
			};
        smart_ptr<Payoff>					payoff_;
        smart_ptr<tree_type>				tree_;
	};


/* A wrapper utility to ease code integration
*/
template<class tree_type, class payoff_type>
class OTS_American_Pricer
{
public:
	typedef tree_type											tree_type;
	typedef typename tree_type::stochastic_process_type			stochastic_process_type;
	typedef template_time_function_utility<payoff_type>				payoff_type;
	typedef OS_American_Pricer_CN<tree_type,payoff_type>		one_step_algorithm_type;
	typedef Cauchy_American_Pricer<tree_type,payoff_type,one_step_algorithm_type>
																cauchy_type;
	OTS_American_Pricer(smart_ptr<stochastic_process_type>& process, 
						smart_ptr<payoff_type>& payoff,
						unsigned long Steps) :
		process_(process),
		payoff_(payoff),
		tree_(new tree_type(process_)),
		cauchy_(tree_,payoff_)
	{
		tree_->set_steps(Steps);
	};
	virtual double operator()() 
	{
			cauchy_.initial_conditions();
			OS_double maturity = payoff_->get_maturity();
			OS_double delta_t = maturity / tree_->get_steps();
			OS_double time_begin = 0;
			OS_double time_end = maturity - delta_t;
			OS_double compound = process_->compound(maturity);
			cauchy_.reset_time();
			cauchy_.insert_time(time_begin);
			cauchy_.insert_time(time_end);

		
			for (;time_begin < time_end; time_begin = ++ cauchy_ ) {};

			OS_double price = compound*jmmath::moment1(cauchy_.get_solution());
			price = std::max(price, (*payoff_)(process_->initial_conditions(),time_begin) );
			return price;
	};

	smart_ptr<stochastic_process_type>				process_;
	smart_ptr<tree_type>							tree_;
	smart_ptr<payoff_type>							payoff_;
	cauchy_type										cauchy_;
};

/* An integration example
*/
template <class payoff_type, class tree_type = OS_tree<geometric_process<double>> >
class OS_Geometric_American_Pricer : public OTS_American_Pricer<tree_type, payoff_type>
{
public:
	typedef OTS_American_Pricer Base;
	OS_Geometric_American_Pricer(double Spot,
									  double Rate,
									  unsigned long Steps,
									  double Expiry,
									  double Volatility,
									  double Strike) :
		Base (
			smart_ptr<stochastic_process_type>(new stochastic_process_type(Rate, Volatility, Spot)),
			smart_ptr<payoff_type>(new payoff_type(Strike,Expiry)),
			Steps
			)
		{};
};



#endif